Libraries:
# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)
library(ggeffects)
library(performance)
Set working directory
knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')
Source sapling data:
source("Scripts/SA_Import.R")
Pivot wider to create dataframe where each row is for one plot and has total details for each species group
sa_merge2 <- sapling_merge %>%
pivot_wider(names_from = Species_Groups, values_from = Total_Tally)
Import time since data and add it to the PIRI dataset
time_since <- read_csv("CleanData/Treat_Year_Data.csv")
sa_merge3 <- merge(sa_merge2, time_since, by = 'Site')
#log transform time from treatment data
sa_merge3$l.TFT <- log(sa_merge3$Time_from_Treat)
Run the ‘Add_BA’ script and merge with dataset:
source("Scripts/Add_BA.R")
# merge with sapling dataset -------------------
sa_merge4 <- merge(sa_merge3, prism_BA, by = 'Plot_No')
Run ‘Ground_Data.R’ script and add it to sapling dataset:
source("Scripts/Ground_Data.R")
# merge with sapling dataset -------------------
sa_merge5 <- merge(sa_merge4, ground3, by = 'Plot_No')
Source and add veg data:
source("Scripts/Veg_Data.R")
# merge with sapling dataset -------------------
sa.all <- merge(sa_merge5, veg3, by = "Plot_No")
rm(sa_merge3,
sa_merge4,
sa_merge5)
Data will be standardized for comparisons at multiple scales
Create a dataframe keeping sapling at 25m2 observation levels & log transforming them
sa.all3 <- sa.all
sa.all3$piri.ba_s <- scale(sa.all3$PIRI.BA_HA)
sa.all3$avgld_s <- scale(sa.all3$avgLD)
sa.all3$vegt_s <- scale(sa.all3$Veg_Total)
sa.all3$min_s <- scale(sa.all3$Mineral_Soil)
sa.all3$so_s <- scale(sa.all3$Shrub_Oak)
sa.all3$other_s <- scale(sa.all3$Other)
sa.all3$ba.ha_s <- scale(sa.all3$BA_HA)
sa.all3$piri_s <- scale(sa.all3$PIRI)
sa.all3 <- sa.all3 %>%
select(Treat_Type, Region, Site, Plot_No, PIRI, piri_s, Shrub_Oak, so_s, Other, other_s, Time_from_Treat, l.TFT, BA_HA, ba.ha_s, PIRI.BA_HA, piri.ba_s, Mineral_Soil, min_s, Litter_Duff, avgLD, avgld_s, Veg_Total, vegt_s) %>%
arrange(Treat_Type)
Looking at pairwise with sa.all3 dataset (sapling counts at 25m2 observation level)
#not log transformed (25m2)
sa3.num <- sa.all3 %>%
select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)
ggpairs(sa3.num)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggpairs(sa3.num, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
rm(sa3.num)
No strong correlations
major personal breakthrough for dummies. For zero inflation, lets see if any TT or region is lacking alltogether in PIRI or SO saplings
sa.all4 <- sa.all3 %>%
group_by(Region)
tapply(sa.all4$PIRI, sa.all4$Region, summary)
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.02459 0.00000 2.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2013 0.0000 9.0000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.07273 0.00000 1.00000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.09574 0.00000 3.00000
# no PIRI in LI,
tapply(sa.all4$PIRI, sa.all4$Treat_Type, summary)
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2564 0.0000 9.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.04902 0.00000 2.00000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.05357 0.00000 2.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.05882 0.00000 2.00000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01149 0.00000 1.00000
I have previously done model elimination with zero inflated poisson distro & with model w.o TT. Neither are better than zi nbinom2 models. zi ~1 has distribution issues, ~Region does not
Model has much better fit and zero inflation issues are resolved when treatment type variable is included. Time to do elimination again but with TT; all models fail with nbinom1 distro, but works with nbinom2 (these notes are from when I was using a zi poisson distro, using nbinom TT no longer seems to be important to model fit, but we will keep for the overarching story)
summary(sa.all3)
## Treat_Type Region Site Plot_No
## Length:498 Length:498 Length:498 Min. : 70.0
## Class :character Class :character Class :character 1st Qu.: 359.5
## Mode :character Mode :character Mode :character Median : 614.0
## Mean : 613.5
## 3rd Qu.: 887.5
## Max. :1143.0
## PIRI piri_s.V1 Shrub_Oak so_s.V1
## Min. :0.00000 Min. :-0.157452 Min. :0.0000 Min. :-0.273816
## 1st Qu.:0.00000 1st Qu.:-0.157452 1st Qu.:0.0000 1st Qu.:-0.273816
## Median :0.00000 Median :-0.157452 Median :0.0000 Median :-0.273816
## Mean :0.09438 Mean : 0.000000 Mean :0.2691 Mean : 0.000000
## 3rd Qu.:0.00000 3rd Qu.:-0.157452 3rd Qu.:0.0000 3rd Qu.:-0.273816
## Max. :9.00000 Max. :14.857482 Max. :7.0000 Max. : 6.849481
## Other other_s.V1 Time_from_Treat l.TFT
## Min. :0.0000 Min. :-0.345712 Min. : 1.0 Min. :0.000
## 1st Qu.:0.0000 1st Qu.:-0.345712 1st Qu.: 3.0 1st Qu.:1.099
## Median :0.0000 Median :-0.345712 Median : 5.0 Median :1.609
## Mean :0.3514 Mean : 0.000000 Mean :10.9 Mean :1.804
## 3rd Qu.:0.0000 3rd Qu.:-0.345712 3rd Qu.: 8.0 3rd Qu.:2.079
## Max. :9.0000 Max. : 8.508476 Max. :33.0 Max. :3.497
## BA_HA ba.ha_s.V1 PIRI.BA_HA piri.ba_s.V1
## Min. : 0.00 Min. :-1.4237625 Min. : 0.0 Min. :-1.285691
## 1st Qu.: 7.00 1st Qu.:-0.8185636 1st Qu.: 5.0 1st Qu.:-0.775675
## Median :16.00 Median :-0.0404508 Median :11.0 Median :-0.163656
## Mean :16.47 Mean : 0.0000000 Mean :12.6 Mean : 0.000000
## 3rd Qu.:23.00 3rd Qu.: 0.5647481 3rd Qu.:18.0 3rd Qu.: 0.550367
## Max. :51.00 Max. : 2.9855436 Max. :44.0 Max. : 3.202451
## Mineral_Soil min_s.V1 Litter_Duff avgLD
## Min. : 0.0000 Min. :-0.197640 Min. : 0.00 Min. : 0.000
## 1st Qu.: 0.0000 1st Qu.:-0.197640 1st Qu.:90.50 1st Qu.: 2.250
## Median : 0.0000 Median :-0.197640 Median :90.50 Median : 3.875
## Mean : 0.7199 Mean : 0.000000 Mean :88.78 Mean : 4.179
## 3rd Qu.: 0.0000 3rd Qu.:-0.197640 3rd Qu.:90.50 3rd Qu.: 5.500
## Max. :50.5000 Max. :13.666913 Max. :90.50 Max. :13.250
## avgld_s.V1 Veg_Total vegt_s.V1
## Min. :-1.652607 Min. : 4.00 Min. :-1.722482
## 1st Qu.:-0.762793 1st Qu.: 37.50 1st Qu.:-0.738485
## Median :-0.120151 Median : 58.00 Median :-0.136337
## Mean : 0.000000 Mean : 62.64 Mean : 0.000000
## 3rd Qu.: 0.522492 3rd Qu.: 83.38 3rd Qu.: 0.609004
## Max. : 3.587405 Max. :202.00 Max. : 4.093383
library(caret)
## Loading required package: lattice
##
## Attaching package: 'caret'
## The following object is masked from 'package:purrr':
##
## lift
findLinearCombos(sa.all3[,4:23]) #shows no linear combos. Don't get how the rank deficiency is happening
## $linearCombos
## $linearCombos[[1]]
## [1] 5 2 3 4
##
## $linearCombos[[2]]
## [1] 7 2 3 6
##
## $linearCombos[[3]]
## [1] 11 2 3 10
##
## $linearCombos[[4]]
## [1] 13 2 3 12
##
## $linearCombos[[5]]
## [1] 15 2 3 14
##
## $linearCombos[[6]]
## [1] 18 2 3 17
##
## $linearCombos[[7]]
## [1] 20 2 3 19
##
##
## $remove
## [1] 5 7 11 13 15 18 20
tapply(sa.all3$Other, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3443 0.0000 6.0000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.411 0.000 6.000
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2143 0.0000 6.0000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.6727 0.0000 9.0000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3511 0.0000 4.0000
tapply(sa.all3$Shrub_Oak, sa.all3$Region, summary) #none in LI
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 1.964 3.000 7.000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2447 0.0000 6.0000
tapply(sa.all3$BA_HA, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.98 23.00 51.00
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 11.00 21.00 19.81 28.00 39.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 14.00 15.32 21.00 44.00
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 13.50 23.00 22.13 31.00 44.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.37 21.00 48.00
tapply(sa.all3$PIRI.BA_HA, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 2 9 10 16 41
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.00 11.00 16.00 17.41 23.00 37.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 9.00 11.44 16.00 34.00
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 8.00 18.00 19.18 30.00 44.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 9.00 10.32 16.00 32.00
tapply(sa.all3$Veg_Total, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 51.50 71.00 70.96 91.88 166.50
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.00 23.00 39.00 44.52 61.00 128.00
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.50 33.75 52.75 56.93 74.38 161.50
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.50 41.25 64.00 62.97 79.75 114.00
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 46.00 62.25 75.09 96.62 202.00
tapply(sa.all3$Mineral_Soil, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.586 0.500 30.500
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.589 0.500 15.500
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.6623 0.0000 50.5000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01818 0.00000 0.50000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2021 0.0000 8.0000
tapply(sa.all3$avgLD, sa.all3$Region, summary) #all present
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.000 1.750 1.977 2.750 7.250
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.75 4.00 5.25 5.61 6.75 12.50
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.250 3.562 5.000 5.347 6.750 13.250
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 2.500 3.500 3.686 4.500 9.750
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 2.812 4.000 4.298 5.438 12.250
tapply(sa.all3$Other, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.8889 1.0000 6.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4118 0.0000 9.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.1429 0.0000 2.0000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.06618 0.00000 2.00000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.1379 0.0000 2.0000
tapply(sa.all3$Shrub_Oak, sa.all3$Treat_Type, summary) #none in mowrx or springrx
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.7094 0.0000 7.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4902 0.0000 6.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
tapply(sa.all3$BA_HA, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 18.00 28.00 26.78 37.00 51.00
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 16.00 17.44 23.00 44.00
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 12.11 18.00 34.00
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 7.000 9.059 14.000 37.000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 10.00 14.00 15.85 21.00 41.00
tapply(sa.all3$PIRI.BA_HA, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 18.00 17.05 25.00 41.00
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 5.00 11.00 14.45 20.25 44.00
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 5.000 9.000 9.804 14.500 30.000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 7.000 7.816 11.000 37.000
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00 9.00 14.00 13.75 18.00 34.00
tapply(sa.all3$Veg_Total, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.00 26.50 46.00 46.85 61.00 115.50
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.5 46.0 68.0 72.7 91.0 202.0
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.00 36.12 55.00 55.21 68.62 132.00
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.00 37.00 66.75 70.47 94.38 196.50
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 16.00 49.50 61.50 64.65 81.00 131.00
tapply(sa.all3$Mineral_Soil, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.008547 0.000000 0.500000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 0.201 0.000 8.000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.652 0.000 50.500
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 0.000 1.199 0.500 30.500
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.9368 0.5000 15.5000
tapply(sa.all3$avgLD, sa.all3$Treat_Type, summary) #all present
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500 3.250 4.750 5.024 6.750 12.500
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.750 3.250 4.375 4.529 5.750 9.750
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.250 3.688 4.775 5.291 6.750 11.000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.500 2.750 2.816 3.812 13.250
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500 1.500 3.500 4.046 5.250 11.500
Best Model AIC: 245.8
sa.m1s <- glmmTMB(PIRI ~ Treat_Type + piri.ba_s + avgld_s + vegt_s + min_s + so_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; false convergence (8). See vignette('troubleshooting'),
## help('diagnose')
#won't converge
#test piri ba
sa.m2s <- glmmTMB(PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m2s) #248
## [1] 247.957
# test ba/ha
sa.m3s <- glmmTMB(PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m3s) #251.7
## [1] 251.6733
lrtest(sa.m2s, sa.m3s) #keep ba
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s +
## ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 19 -104.98
## 2 18 -107.84 -1 5.7163 0.01681 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test SO
sa.m4s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m4s) #248
## [1] 248.2719
lrtest(sa.m4s, sa.m2s) #drop
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + avgld_s + vegt_s + min_s + so_s + other_s +
## ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 18 -106.14
## 2 19 -104.98 1 2.3148 0.1281
# test other
sa.m5s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m5s) #251.1
## [1] 251.0924
lrtest(sa.m4s, sa.m5s) #keep other
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 18 -106.14
## 2 17 -108.55 -1 4.8206 0.02812 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test mineral
sa.m6s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m6s) #247
## [1] 247.0116
lrtest(sa.m4s, sa.m6s) #drop
## Likelihood ratio test
##
## Model 1: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + min_s + other_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + ba.ha_s + avgld_s + vegt_s + other_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 18 -106.14
## 2 17 -106.51 -1 0.7397 0.3897
# test veg
sa.m7s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + avgld_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m7s) #245.8
## [1] 245.8436
# test ld
sa.m8s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m8s) #244.8
## [1] 244.7564
sa.m8s_sr <- simulateResiduals(sa.m8s, n = 1000, plot = TRUE) #passes
testResiduals(sa.m8s_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.029239, p-value = 0.7881
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.4242, p-value = 0.826
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.029239, p-value = 0.7881
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.4242, p-value = 0.826
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00072289156626506 )
## 0
# passes
emmeans(sa.m8s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response') # the model is rank deficient
## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.0272 NaN Inf NaN NaN
## FallRx 0.0235 NaN Inf NaN NaN
## Harvest 0.1676 NaN Inf NaN NaN
## MowRx 0.1509 NaN Inf NaN NaN
## SpringRx 0.0193 NaN Inf NaN NaN
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 1.157 NaN Inf 1 NaN NaN
## Control / Harvest 0.162 NaN Inf 1 NaN NaN
## Control / MowRx 0.180 NaN Inf 1 NaN NaN
## Control / SpringRx 1.409 NaN Inf 1 NaN NaN
## FallRx / Harvest 0.140 NaN Inf 1 NaN NaN
## FallRx / MowRx 0.156 NaN Inf 1 NaN NaN
## FallRx / SpringRx 1.217 NaN Inf 1 NaN NaN
## Harvest / MowRx 1.111 NaN Inf 1 NaN NaN
## Harvest / SpringRx 8.678 NaN Inf 1 NaN NaN
## MowRx / SpringRx 7.812 NaN Inf 1 NaN NaN
##
## P value adjustment: tukey method for comparing a family of 2 estimates
## Tests are performed on the log scale
sa.m9s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m9s) #245.9 with veg; 245.8 with ld
## [1] 245.8436
summary(sa.m9s)
## Family: nbinom2 ( log )
## Formula:
## PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## Zero inflation: ~Region
## Data: sa.all3
##
## AIC BIC logLik deviance df.resid
## 245.8 313.2 -106.9 213.8 482
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.00000003136 0.0001771
## Site (Intercept) 1.00553141354 1.0027619
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for nbinom2 family (): 0.332
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.4223 0.8811 -6.154 0.000000000755 ***
## Treat_TypeFallRx -0.2051 0.9401 -0.218 0.8273
## Treat_TypeHarvest 1.8988 1.3009 1.460 0.1444
## Treat_TypeMowRx 1.6541 1.0395 1.591 0.1116
## Treat_TypeSpringRx -0.2551 1.4500 -0.176 0.8603
## ba.ha_s 0.8721 0.3623 2.407 0.0161 *
## other_s -1.1511 0.6811 -1.690 0.0910 .
## avgld_s -0.2776 0.2909 -0.954 0.3398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.521 1.019 1.492 0.136
## RegionLI 22.406 53458.279 0.000 1.000
## RegionMA -20.529 22550.145 -0.001 0.999
## RegionME -2.795 3.323 -0.841 0.400
## RegionNH -1.687 1.457 -1.158 0.247
Model without avgld is rank deficient
sa.m9s <- glmmTMB(PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom2)
AIC(sa.m9s) #245.9 with veg; 245.8 with ld
## [1] 245.8436
summary(sa.m9s)
## Family: nbinom2 ( log )
## Formula:
## PIRI ~ Treat_Type + ba.ha_s + other_s + avgld_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## Zero inflation: ~Region
## Data: sa.all3
##
## AIC BIC logLik deviance df.resid
## 245.8 313.2 -106.9 213.8 482
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.00000003136 0.0001771
## Site (Intercept) 1.00553141354 1.0027619
## Number of obs: 498, groups: Plot_No:Site, 498; Site, 47
##
## Dispersion parameter for nbinom2 family (): 0.332
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.4223 0.8811 -6.154 0.000000000755 ***
## Treat_TypeFallRx -0.2051 0.9401 -0.218 0.8273
## Treat_TypeHarvest 1.8988 1.3009 1.460 0.1444
## Treat_TypeMowRx 1.6541 1.0395 1.591 0.1116
## Treat_TypeSpringRx -0.2551 1.4500 -0.176 0.8603
## ba.ha_s 0.8721 0.3623 2.407 0.0161 *
## other_s -1.1511 0.6811 -1.690 0.0910 .
## avgld_s -0.2776 0.2909 -0.954 0.3398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.521 1.019 1.492 0.136
## RegionLI 22.406 53458.279 0.000 1.000
## RegionMA -20.529 22550.145 -0.001 0.999
## RegionME -2.795 3.323 -0.841 0.400
## RegionNH -1.687 1.457 -1.158 0.247
sa.m9s_sr <- simulateResiduals(sa.m9s, n = 1000, plot = TRUE) #passes
testResiduals(sa.m9s_sr)
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.030633, p-value = 0.7383
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.30212, p-value = 0.862
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000883534136546185 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.030633, p-value = 0.7383
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.30212, p-value = 0.862
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000883534136546185 )
## 0
# passes
emmeans(sa.m9s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.0268 0.0236 Inf 0.00477 0.151
## FallRx 0.0219 0.0206 Inf 0.00345 0.139
## Harvest 0.1792 0.1821 Inf 0.02444 1.314
## MowRx 0.1403 0.1040 Inf 0.03279 0.600
## SpringRx 0.0208 0.0280 Inf 0.00148 0.292
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 1.228 1.154 Inf 1 0.218 0.9995
## Control / Harvest 0.150 0.195 Inf 1 -1.460 0.5889
## Control / MowRx 0.191 0.199 Inf 1 -1.591 0.5030
## Control / SpringRx 1.291 1.871 Inf 1 0.176 0.9998
## FallRx / Harvest 0.122 0.166 Inf 1 -1.542 0.5352
## FallRx / MowRx 0.156 0.170 Inf 1 -1.704 0.4314
## FallRx / SpringRx 1.051 1.577 Inf 1 0.033 1.0000
## Harvest / MowRx 1.277 1.562 Inf 1 0.200 0.9996
## Harvest / SpringRx 8.619 14.264 Inf 1 1.302 0.6903
## MowRx / SpringRx 6.748 10.028 Inf 1 1.285 0.7008
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale
None significantly different from eachother
sa.piri1 <- ggpredict(sa.m9s, terms = c("avgld_s", "Treat_Type"))
#control and fallrx are almost completely identical - if I wanted to nudge/jitter
# sa.piri2 <- sa.piri1 %>%
# filter(group == "FallRx")
#
# sa.piri2$predicted <- sa.piri2$predicted+0.01
#
# sa.piri3 <- sa.piri1 %>%
# filter(group != "FallRx")
#
# sa.piri3 <- rbind(sa.piri3, sa.piri2)
ggplot(sa.piri1, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
labs(x = "Average leaf litter depth \n (standardized)",
y = NULL)
ggsave(filename = "SA_PIRI_avgld.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)
sa.piri4 <- ggpredict(sa.m9s, terms = c("ba.ha_s", "Treat_Type"))
#again control and fall rx overlap - if I wanted to nudge/jitter
# sa.piri5 <- sa.piri4 %>%
# filter(group == "FallRx")
#
# sa.piri5$predicted <- sa.piri5$predicted+0.01
#
# sa.piri6 <- sa.piri4 %>%
# filter(group != "FallRx")
#
# sa.piri6 <- rbind(sa.piri6, sa.piri5)
ggplot(sa.piri4, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
labs(x = "Average basal area per hectare \n (standardized)",
y = "Pitch pine stems \n (adjusted for time)")
ggsave(filename = "SA_PIRI_BA.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)
sa.piri5 <- ggpredict(sa.m9s, terms = c("other_s", "Treat_Type"))
ggplot(sa.piri5, aes(x, predicted, color = group))+
geom_line(linewidth = 1, show.legend = FALSE) +
scale_color_manual(values = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))+
theme_classic()+
theme(panel.background = element_blank()) +
labs(colour = "Treatment Type")+
theme(axis.text = element_text(size = 18), axis.title = element_text(size = 24))+
labs(x = "Counts of other saplings \n (standardized)",
y = NULL)
ggsave(filename = "SA_PIRI_other.tiff", path = "Plots/PIRI_GLMM", width = 6, height = 4.2, device = "tiff", dpi = 700)
Best Model AIC: 323.8
Start from the beginning knowing that both Region and Treat Type have 0 areas
tapply(sa.all4$Shrub_Oak, sa.all4$Region, summary)
## $ALB
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01639 0.00000 1.00000
##
## $LI
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $MA
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.000000 0.000000 0.006494 0.000000 1.000000
##
## $ME
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.000 1.000 1.964 3.000 7.000
##
## $NH
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.2447 0.0000 6.0000
# no SO in LI,
tapply(sa.all4$Shrub_Oak, sa.all4$Treat_Type, summary)
## $Control
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.7094 0.0000 7.0000
##
## $FallRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4902 0.0000 6.0000
##
## $Harvest
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.00000 0.00000 0.00000 0.01786 0.00000 1.00000
##
## $MowRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
##
## $SpringRx
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0
# no SO in MowRX or SpringRx
rm(sa.all4)
sa.all_noSMRX <- sa.all3 %>%
filter(Treat_Type != "SpringRx") %>%
filter(Treat_Type != "MowRx")
There is just not enough information in exposed mineral soil (so many 0s) for it to be useful. It messes up models, making them rank deficient. So just dropping that variable
so.sa1s <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + avgld_s + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
# won't converge
#had to remove piri ba for model to converge
so.sa2s <- glmmTMB(Shrub_Oak ~ Treat_Type + avgld_s + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa2s) #333
## [1] 333.03
#avgld
so.sa3s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa3s) #331
## [1] 331.278
lrtest(so.sa2s, so.sa3s) #0.6
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + avgld_s + vegt_s + piri_s + other_s +
## ba.ha_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 16 -150.51
## 2 15 -150.64 -1 0.248 0.6185
#piri
so.sa4s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa4s) #329
## [1] 329.3011
lrtest(so.sa3s, so.sa4s) #0.9
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + piri_s + other_s + ba.ha_s +
## offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 15 -150.64
## 2 14 -150.65 -1 0.0231 0.8793
#other
so.sa5s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa5s) #327.3
## [1] 327.3061
lrtest(so.sa4s, so.sa5s) #0.9 - drop other
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + other_s + ba.ha_s + offset(l.TFT) +
## (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 14 -150.65
## 2 13 -150.65 -1 0.005 0.9436
#test BA
so.sa6s <- glmmTMB(Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC (so.sa6s) #325.3
## [1] 325.3149
lrtest(so.sa6s, so.sa5s) #0.9
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + ba.ha_s + offset(l.TFT) + (1 |
## Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 12 -150.66
## 2 13 -150.65 1 0.0088 0.9253
# test veg
so.sa7s <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa7s) #323.8
## [1] 323.7947
lrtest(so.sa7s, so.sa6s) #0.5
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + vegt_s + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -150.90
## 2 12 -150.66 1 0.4798 0.4885
# add piri ba back in
so.sa8s <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
AIC(so.sa8s) #325
## [1] 325.0044
lrtest(so.sa7s, so.sa8s) #0.4
## Likelihood ratio test
##
## Model 1: Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1 | Site/Plot_No)
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 11 -150.9
## 2 12 -150.5 1 0.7902 0.374
rm(so.sa1s,
so.sa2s,
so.sa3s,
so.sa4s,
so.sa5s,
so.sa6s,
so.sa8s)
so.sa7s <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all_noSMRX,
ziformula = ~Region,
family = nbinom1)
summary(so.sa7s) #323.8
## Family: nbinom1 ( log )
## Formula: Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Zero inflation: ~Region
## Data: sa.all_noSMRX
##
## AIC BIC logLik deviance df.resid
## 323.8 363.6 -150.9 301.8 264
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Plot_No:Site (Intercept) 0.000000007557 0.00008693
## Site (Intercept) 0.109978958497 0.33163076
## Number of obs: 275, groups: Plot_No:Site, 275; Site, 26
##
## Dispersion parameter for nbinom1 family (): 0.773
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6952 0.2583 -10.435 <0.0000000000000002 ***
## Treat_TypeFallRx 0.8099 0.3689 2.195 0.0281 *
## Treat_TypeHarvest 1.5834 1.2203 1.298 0.1944
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9034 0.7926 2.402 0.0163 *
## RegionLI 30.0378 1411573.3749 0.000 1.0000
## RegionMA 2.0029 1.3083 1.531 0.1258
## RegionME -31.2661 873734.0264 0.000 1.0000
## RegionNH -1.0435 0.8939 -1.167 0.2431
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
so.sa7s_sr <- simulateResiduals(so.sa7s, n = 1000, plot = T) #passes
testResiduals(so.sa7s_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.032251, p-value = 0.9372
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9179, p-value = 0.902
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 275, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.005545455
## sample estimates:
## outlier frequency (expected: 0.00112727272727273 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.032251, p-value = 0.9372
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.9179, p-value = 0.902
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 275, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.005545455
## sample estimates:
## outlier frequency (expected: 0.00112727272727273 )
## 0
emmeans(so.sa7s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.722 0.187 Inf 0.435 1.20
## FallRx 1.623 0.590 Inf 0.796 3.31
## Harvest 3.518 4.244 Inf 0.331 37.42
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio p.value
## Control / FallRx 0.445 0.164 Inf 1 -2.195 0.0719
## Control / Harvest 0.205 0.250 Inf 1 -1.298 0.3965
## FallRx / Harvest 0.461 0.570 Inf 1 -0.626 0.8056
##
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scale
Control vs FallRx close at 0.07
Need to figure out better use of emmeans (potentially)
so.sa1.plot <- emmeans(so.sa7s, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
plot <- plot(so.sa1.plot, horizontal=FALSE, comparisons=TRUE, color = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15")) + theme_bw()
plot+
theme_classic()+
theme(panel.background = element_blank()) +
theme(axis.text = element_text(size = 14), axis.title = element_text(size = 20), legend.position = "right", legend.text = element_text(size = 14), legend.title = element_text(size = 18)) +
labs(x = "Shrub Oak Stems \n (corrected for time from treatment)",
z = "Treatment Type",
color = "Treatment Type")
More in harvest - does fire kill large scrub oak saplings? They are entirely missing from springrx and mowrx treatments
There are no shrub oak sapling in springrx and mowrx treatments, so the above model is better - but just to have it.
so.all1 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri_s + other_s + min_s + avgld_s+ piri.ba_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
## Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
## NA/NaN function evaluation
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; function evaluation limit reached without convergence (9). See
## vignette('troubleshooting'), help('diagnose')
#won't converge
#drop min & piri ba
so.all2 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri_s + other_s + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all2) #335
## [1] 335.3884
so.all3 <- glmmTMB(Shrub_Oak ~ Treat_Type + other_s + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all3) #333.5
## [1] 333.478
so.all4 <- glmmTMB(Shrub_Oak ~ Treat_Type + avgld_s + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all4) #331.6
## [1] 331.5949
so.all5 <- glmmTMB(Shrub_Oak ~ Treat_Type + ba.ha_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all5) #329.7
## [1] 329.725
so.all6 <- glmmTMB(Shrub_Oak ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all6) #327.8
## [1] 327.7947
so.all7 <- glmmTMB(Shrub_Oak ~ Treat_Type + piri.ba_s + offset(l.TFT) + (1|Site/Plot_No),
data = sa.all3,
ziformula = ~Region,
family = nbinom1)
AIC(so.all7) #329
## [1] 329.0044
rm(so.all1,
so.all2,
so.all3,
so.all4,
so.all5,
so.all7)
so.all6_sr <- simulateResiduals(so.all6, n = 1000, plot = T) #passes
testResiduals(so.all6_sr) #passes
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.026629, p-value = 0.8718
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.88884, p-value = 0.844
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000522088353413655 )
## 0
## $uniformity
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.026629, p-value = 0.8718
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.88884, p-value = 0.844
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000522088353413655 )
## 0
emmeans(so.all6, specs = pairwise ~ Treat_Type, adjust = 'Tukey', type = 'response')
## $emmeans
## Treat_Type response SE df asymp.LCL asymp.UCL
## Control 0.410 0.1059479 Inf 0.247 1
## FallRx 0.922 0.3350685 Inf 0.452 2
## Harvest 1.998 2.4106105 Inf 0.188 21
## MowRx 0.000 0.0000001 Inf 0.000 Inf
## SpringRx 0.000 0.0000044 Inf 0.000 Inf
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scale
##
## $contrasts
## contrast ratio SE df null z.ratio
## Control / FallRx 0 0 Inf 1 -2.196
## Control / Harvest 0 0 Inf 1 -1.298
## Control / MowRx 10634534517938 16507512852342061056 Inf 1 0.000
## Control / SpringRx 10155053986 1101249307522083 Inf 1 0.000
## FallRx / Harvest 0 1 Inf 1 -0.626
## FallRx / MowRx 23902818311903 37103277047566393344 Inf 1 0.000
## FallRx / SpringRx 22825109079 2475234065907171 Inf 1 0.000
## Harvest / MowRx 51809117272110 80420978256787030016 Inf 1 0.000
## Harvest / SpringRx 49473193394 5365044838234219 Inf 1 0.000
## MowRx / SpringRx 0 1486 Inf 1 0.000
## p.value
## 0.1813
## 0.6927
## 1.0000
## 1.0000
## 0.9709
## 1.0000
## 1.0000
## 1.0000
## 1.0000
## 1.0000
##
## P value adjustment: tukey method for comparing a family of 5 estimates
## Tests are performed on the log scale